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The  present  paper  describes  an  algorithm  for  the  rapid  evaluation  of  the  potential  and  force  fields 
in  systems  involving  large  numbers  of  particles  whose  interactions  are  described  by  Coulomb’s  law. 
Unlike  previously  published  schemes,  the  algorithm  of  this  paper  has  an  asymptotic  CPU  time 
estimate  of  0(7V),  where  N  is  the  number  of  particles  in  the  simulation,  and  does  not  depend  for 
Its  efficient  performance  on  the  statbtics  of  the  distribution.  The  numerical  examples  we  present 
indicate  that  it  should  be  an  algorithm  of  choice  in  many  situations  of  practical  interest. 
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1.  Introduction 

The  evaluation  of  Coulombic  and  gravitational  interactions  in  large-scale  ensembles  of  particles 
is  an  integral  part  of  the  numerical  simulation  of  a  large  number  of  physical  processes.  Ty’pical 
examples  include  celestial  mechanics,  plasma  simulations,  the  vortex  method  in  fluid  dynamics, 
and  the  solution  of  the  Laplace  equation  via  potential  theory  (see  (l],  [2],  [3],  [8],  [10]).  In  such 
cases,  the  potential  has  the  form 

^  “  ^external  ^locat  "b  ^/ar,  (l) 

where  ^loeat  is  a  rapidly  decaying  function  of  distance  (such  as  the  Van  der  Waals  potential 
in  chemical  physics),  ^eiUmaf  is  a-fimction  which  is  independent  of  the  number  and  relative 
positions  of  the  particles  (such  as  an  external  gravitational  field),  and  9 far  is  Coulombic  or 
gravitational. 

In  the  numerical  evaluation  of  fields  of  the  form  (1),  the  cost  of  computing  the  terms 
^external  nnd  ^toeat  is  of  the  order  O(A'),  where  N  is  the  number  of  particles  in  the  ensemble. 
Indeed,  ^txtemai  is  evaluated  separately  for  each  particle,  and  9ioeai  decays  rapidly,  involving 
the  interactions  of  each  particle  with  a  small  number  of  nearest  neighbors.  Unfortunately,  eval¬ 
uation  of  the  term  4»/ar,  if  done  directly,  requires  order  O(N^)  operations,  since  the  Coulombic 
potential  decays  slowly,  and  the  interactions  between  each  pair  of  particles  have  to  be  taken 
into  account.  In  many  situations,  in  order  to  be  of  physical  interest,  the  simulation  has  to 
involve  thousands  of  particles  (or  more),  making  the  estimate  0(A'*)  excessive  in  some  cases, 
and  prohibitive  in  others. 

Several  different  approaches  have  been  used  to  reduce  the  cost  of  the  Coulombic  part  of 
the  computation.  For  a  detailed  discussion  of  these  algorithms,  we  refer  the  reader  to  [7],  and 
to  the  original  papers  [l],  [2],  [8],  (10).  Here,  we  just  observe  that  each  of  the  algorithms  [l], 
[2],  [7],  [8],  [10]  imposes  strong  requirements  on  the  statistics  of  the  charge  distribution.  In 
particular,  the  methods  of  [1],  [7],  and  [8]  require  that  the  distribution  be  reasonably  uniform 
in  a  square-shaped  region  of  interest,  the  algorithm  of  [10]  assumes  that  the  charges  are  located 
on  a  curve  in  and  the  algorithm  of  [2]  works  fa'rly  well  for  highly  clustered  distributions, 
but  faik  for  uniform  ones. 

In  this  paper,  we  introduce  an  algorithm  for  the  rapid  evaluation  of  the  potential  and  force 
fields  of  large-scale  ensembles  of  particles.  To  evaluate  all  Coulombic  interactions  of  N  particles 
in  R^,  the  algorithm  requires  an  amount  of  work  of  the  order  0(A' ),  and  this  estimate  does  not 
depend  on  the  statistics  of  the  distribution. 

The  procedure  described  here  is  an  adaptive  version  of  the  algorithm  of  [7].  In  the  follow¬ 
ing  section,  we  introduce  the  analytical  apparatus  to  be  used.  Section  3  contains  a  detailed 
description  of  the  algorithm  and  its  complexity  analysis,  and  in  Section  4  we  present  numerical 
experiments  demonstrating  the  actual  performance  of  the  scheme. 

Remark  1.1.  Given  a  collection  of  points  zi,'--,Zn  in  C,  the  Hilbert  matrix  associated  with 
the  points  {r,}  is  defined  as  follows: 

Zj  ”  Zj 
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An  =  0. 


It  immediately  follows  from  Lemma  2.1  and  formula  (3)  that  evaluating  the  fields  of  a  set  of 
charges  of  strengths  qi,"  •  ,qn  located  at  the  points  •••,«„  at  these  same  points  is  equivalent 
to  applying  the  associated  Hilbert  matrix  to  the  vector  (91,-  --.ffn).  Therefore,  the  algorithm 
of  the  present  paper  may  be  viewed  as  an  order  0(n)  procedure  for  applying  an  n  x  n  Hilbert 
matrix  to  an  arbitrary  vector.  Recently,  several  papers  have  been  published  on  this  subject, 
also  referred  as  the  TVummer  problem,  (see  |6],  [4],  (5),  [9]). 


2.  Analytical  tools 

VVe  consider  a  two  dimensional  physical  model  consisting  of  a  set  of  particles  whose  pairwise 
interactions  are  described  by  Coulomb’s  law\  More  precisely,  suppose  that  a  point  charge  of 
unit  strength  is  located  at  the  point  Xo  =  (jfo*yo)  €  R^-  Then,  for  any  x  =  (x.y)  €  \  {xo). 

the  potential  and  electrostatic  field  due  to  this  charge  are  described  by  the  expressions: 

^o(x)  = -/o9(I|x-xo||)  (2) 


,W  = 


(x  -  Xq) 
||x  -  Xo||2 


respectively. 

It  is  well  known  that  the  function  <i>xo  is  harmonic  in  any  region  not  containing  Xo.  and 
that  for  every  harmonic  function  «,  there  exists  an  analytic  function  ix? :  C  -*  C  such  that  v  is 
the  real  part  of  w.  In  particular,  we  have 

^o(x)  =  Re{-log(!’  -  Zq)).  (4) 

In  the  remainder  of  this  paper  we  shall  work  with  analytic  functions  in  C,  making  no  distinction 
between  a  point  (or,  y)  €  R^  and  a  point  i  +  ly  €  C.  Following  standard  practice,  we  w-ill  refer 
to  the  analytic  function  log{z)  as  the  potential  due  to  a  charge.  For  more  complicated  charge 
distributions,  we  w’ill  use  other  analytic  functions  and  we  will  refer  to  them  as  potentials  as 
well. 

Detailed  proofs  of  Lemmas  2.1  -  2.4  and  Theorem  2.1  below  can  be  found  in  [7]. 

Lemma  2.1.  If  u(x,y)  =  Re{w{x,y))  describes  the  potential  field  at  (z,y),  then  the  force  field 
is  given  by 

Vu  =  =  (Re(u>'),-/m(w'))» 

where  w'  is  the  derivative  of  w. 

The  following  theorem  gives  the  expression  for  the  multipole  expansion  of  the  potential  due 
to  a  set  of  charges  and  an  estimate  for  the  remainder  of  this  expansion  after  k  terms. 

Theorem  2.1.  Suppose  that  m  charges  of  strengths  9,-,  1 »  !,•••, m  are  located  at  points 
Zi,  I  =  1,  -  •  •  ,m,  with  |2,I  <  r.  Then,  for  any  with  Is)  >  r,  the  potential  ^(z)  is  given  by 
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where 


=  <?  - /^(i)  +  52  (5) 

*-i  * 


anrf 


Furthermore,  for  any  p>  1, 


where 


C=\^  .  ^  =  £k.l  »  anrf  a= 


We  will  use  a  simple  example  to  demonstrate  how  multipole  expansions  can  be  used  to  speed 
up  calculations  with  potentials.  Suppose  that  X  =  {ii^ara^ •  •  •  »3fm}  and  Y  =  {yi,y2,- •  •  ,yn} 
are  two  finite  sets  of  points  in  C.  We  say  that  the  sets  A'  and  Y  are  well  separated  (Figure  1) 
if  there  exist  two  points  zqi  yo  €  <D  and  a  real  number  r  >  0  such  that 


\xj  -  xo\<  r  for  all  i  =  1,- 
ly;  ~  yo|  <  »■  for  all  i  =  1,- 
|xo  -  yo|  >  3r. 

Suppose  now  that  charges  of  strengths  {?i ,  ? j,  •  •  • ,  9m  }  are  located  at  the  points  {xj ,  xj ,  •  •  • ,  i m  } 
and  that  we  wish  to  evaluate  the  sum 

^Myj)  (9) 

•  si 

for  all  j  =  1,2,'  •  •  ,n.  Clearly,  this  requires  order  n  •  m  work  (evaluating  m  fields  at  n  points). 
Now  suppose  that  we  first  evaluate  the  coefficients  of  a  p-tcrm  multipole  expansion  due  to  the 
charges  {71,92,- ••  ,9m}  about  xq,  using  Theorem  2.1.  This  requires  a  number  of  operations 
proportional  to  m  •  p.  Evaluating  the  resulting  multipole  expansion  at  all  points  yy  requires 
order  n-p  work,  and  the  total  computational  effort  is  of  the  order  0(m  pH- n  p).  Furthermore, 
due  to  (7), 

(10) 

and  in  order  to  obtain  a  relative  precision  e,  p  must  be  of  the  order  — /092(e) •  Once  the  precision 
b  specified,  the  amount  of  computation  b  reduced  to  0(n  +  m),  which  is  significantly  smaller 
than  m  •  n  for  large  m  and  n. 

The  following  three  lemmas  describe  translation  operators  for  multipole  and  power  series 
expansions  in  R^,  and  provide  error  bounds  adlowing  the  manipulation  of  these  expansions  in 
the  manner  required  the  algorithm.  The  first.  Lemma  2.2,  supplies  a  mechanism  for  shifting 
the  center  of  a  multipole  expansion. 


£^*.(y/)  -  Qiog[vi  ~  »o)  -  S  7-  ^  ^ 

1-1  *-i  (y/  ~  ^0) 
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Lemma  2.2.  Suppose  that 


tf>{z)  =  aolog{z  -  zo)  +  u 


is  a  multipole  expansion  of  the  potential  due  to  a  set  of  m  charges  of  strengths  gi,q2,’“ » 9m. 
all  of  which  are  located  inside  the  circle  D  of  radius  R  with  center  at  zq.  Then  for  z  outside 
the  circle  D\  of  radius  {R  +  |2oI)  centered  at  the  origin, 


<f>{z)  =  aolog{z)  + 


where 


with  (j[)  the  binomial  coefficients.  Furthermore,  for  any  p  >  1, 

d>{z)  -  aolog{z)  '  Z  7  ^  (~  |Voi^fi|) 
with  A  defined  in  Theorem  2.1. 

Lemma  2.3  describes  the  conversion  of  a  multipole  expansion  into  a  local  (Taylor)  expansion 
in  a  circular  region  of  analyticity. 

Lemma  2.3.  Suppose  that  m  charges  of  strengths  9i.92.  ’  •  •»9m.  ore  located  inside  the  circle 
D\  with  radius  R  and  center  at  zo,  and  that  |zo|  >  (c+  l)i?  with  c>  1.  Then  the  corresponding 
multipole  expansion  (11)  converges  inside  the  circle  D2  of  radius  R  centered  about  the  origin. 
Inside  D2,  the  potential  due  to  the  charges  is  described  by  a  power  series: 

d>{z)  =  '£brz',  (15) 


where 


bo  =  aolog{~zo)  +  ^ 

kml  ^0 
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Furthermore,  for  any  p  >  maa'(2,  ^),  an  error  bound  for  the  truncated  series  is  given  by 


A 


m  - 

where  A  is  defined  in  Theorem  B.l  and  e  is  the  base  of  natural  logarithms. 

Remark  2.1.  In  Theorem  2.1,  the  charges  {quQ2t“'  tQm}  can  be  replaced  with  dipoles,  or 
with  finite  linear  combinations  of  multipoles  of  the  form 

"<>• 

ao/o<7(2)  + (19) 

In  this  case,  the  form  of  the  expressions  (5)-(8)  is  unchanged.  However,  the  coefficients  Q,  A, 
o,  {a*},  k  =  1,2,  •  •  • ,  now  depend  on  oOfOi,’'  •,£»«,  and  can  be  easily  determined  by  repeated 
differentiation  of  (6)  with  respect  to  2,-,  i  =  1,2,‘ •  •,m. 

Remark  2.2.  Ifin  Lemma  2.3,  the  field  <l>{z)  is  generated  by  a  single  charge  located  at  zq, 
then  the  only  non-zero  term  in  the  expansion  (11)  is  oo,  and  oq  =  9i.  Similarly,  if  the  field 
^(2)  is  generated  by  a  single  dipole  located  at  20,  then  the  only  non-zero  term  in  the  expansion 
(11)  is  ai,  and  ai  =  91. 

Lemma  2.4  provides  a  formula  for  shifting  the  center  of  a  local  expansion.  The  expression 
(20)  below  is  an  exact  one,  and  no  error  bound  is  needed. 

Lemma  2.4.  For  any  complex  20,2  and  {a*},  k  =  1,2, 

t  -  «)‘  =  L  ft  f  ? )  *'■  <“) 

*«=0  /«0  \kml  \  /  ) 

Remark  2.3.  One  of  the  advantages  of  using  expansions  of  the  forms  (5)  and  (15)  for  repre¬ 
senting  potential  ^elds  is  the  fact  that  these  expansions  possess  simple  analytical  derivatives. 
This  permits  the  force  fields  to  be  obtained  from  the  potentials  by  Lemma  2.1,  without  the  use 
of  numerical  differentiation  and  the  attendant  loss  of  accuracy. 

3.  Ths  adaptive  multipole  algorithm 

3.1  General  strategy 

In  this  section,  we  describe  an  adaptive  algorithm  for  the  rapid  evaluation  of  the  potential 
and  electrostatic  fields  due  to  arbitrary  distributions  of  charges  and/or  dipoles.  The  main 
strateg>’  b  similar  to  that  described  in  [7].  It  consists  in  clustering  particles  at  dijfferent  spatial 
lengths  and  using  multipole  expansions  to  evaluate  the  interactions  between  clusters  which  are 
sufficiently  far  away  from  each  other.  The  interactions  between  particles  which  are  nearby  are 
computed  directly. 

To  be  more  specific,  consider  the  domain  depicted  in  Figure  2.  N  charges  are  arbitrarily 
dbtributed  in  R’,  and,  without  loss  of  generality,  we  can  assume  that  all  of  them  are  located 
inside  a  square  with  sides  of  length  one,  centered  about  the  origin  of  the  coordinate  system. 
This  square  will  be  referred  to  as  the  computational  box. 


^6,  2'  <  l)  +  c^)  /l^' 

M  '  c(c-l)  Vc/ 
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Given  a  machine  precision  £,  we  set  the  number  of  terms  in  all  expansions  to  p  fa  |/op2(c)|, 
and  specify  that  no  interactions  be  evaluated  via  multipole  expansions  for  groups  of  particles 
that  are  not  well-separated.  This  is  precbely  the  condition  needed  for  the  error  bounds  (7), 
(14),  (18)  to  apply  with  c  =  2.  In  order  to  impose  such  a  condition,  we  introduce  a  hierarchy  of 
meshes  which  refine  the  computational  box  into  smaller  and  smaller  regions  (Figure  3).  Mesh 
level  0  corresponds  to  the  entire  computational  box,  while  mesh  level  /  -b  1  is  obtained  from 
mesh  level  I  by  subdividing  each  region  into  four  equal  parts.  A  tree  structure  is  imposed  on 
this  mesh  hierarchy,  so  that  if  b  is  a  fixed  box  at  level  I,  the  four  boxes  at  level  I  H- 1  obtained 
by  subdividing  b  are  considered  its  children.  The  four  children  of  the  same  box  will  be  referred 
to  as  brothers. 

However,  unlike  the  algorithm  of  (7),  we  do  not  use  the  same  number  of  levels  for  all  parts 
of  the  computational  box.  Instead,  some  integer  5  >  0  is  fixed,  and  at  every  level  of  refinement 
we  subdivide  only  those  boxes  that  contain  more  than  s  charges.  Generally,  this  results  in  a 
large  number  of  empty  boxes  at  finer  levels  of  the  procedure.  At  every  level  of  refinement,  a 
table  of  non-empty  boxes  is  maintained,  so  that  once  an  empty  box  is  encountered,  its  existence 
b  immediately  forgotten  and  it  is  completely  ignored  by  the  subsequent  process. 

Observation  3.1.  It  should  be  noted  that  for  a  fixed  machine  precision  t,  only  certain  classes 
of  particle  distributions  can  be  modeled,  independently  of  the  algorithm  used.  Namely,  suppose 
that  two  charges  ci,  C2  in  a  distribution  have  positions  ii,  12  and  that  |lxi  -  ar2ll  <  5  •  Iki  +  *’2ll* 
Obviously,  under  these  conditions  the  particles  cj,  C2  can  not  be  discerned,  and  no  meaningful 
simulation  b  possible.  Since  the  smallest  discernible  distance  between  tw’o  particles  depends 
on  the  position  of  these  particles  in  the  computational  cell,  such  position-dependent  condition 
can  not  be  imposed  a  priori.  In  order  to  make  the  simulation  possible,  we  will  simply  require 
that  r„„n  >  c,  where  r^in  is  the  smallest  distance  between  any  two  particles  in  the  simulation, 
and  €  is  the  machine  precision.  Therefore,  the  maximum  number  of  ancestors  for  any  box  in 
the  computational  cell  is  p=  ||Iop2(«)l  • 

3.2  Notation 

In  this  subsection,  we  introduce  several  definitions  to  be  used  in  the  description  of  the  algorithm 
below. 

For  any  subset  A  of  the  computational  box,  T{A)  will  denote  the  set  of  particles  that  are 
contained  in  A. 

:  i^  tht  set  of  non-empty  boxes  at  level  t.  Bo  consists  of  only  the  computational  box  itself. 
vVe  will  denote  by  nlev  the  highest  level  of  refinement  at  any  point. 

If  a  box  contains  more  than  9  particles,  it  b  called  a  parent  box.  Otherwise,  the  box  is  said  to 
be  ehildleas. 

A  child  box  b  a  non-empty  box  resulting  from  the  division  into  four  of  a  parent  box. 
Colleaguea  are  adjacent  boxes  of  the  same  size  (at  the  same  level).  A  given  box  has  at  most  8 
colleagues  (Figure  4). 

With  each  box  6  at  level  I  we  will  associate  five  lists  of  other  boxes,  determined  by  their 
positions  with  respect  to  b.  Following  are  the  definitions  of  these  lists  (Figure  5). 
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List  i  of  a  box  6  will  be  denoted  by  17»;  it  is  empty  if  6  is  a  parent  box.  If  6  is  childless,  Ui 
consists  of  b  and  all  childless  boxes  adjacent  to  b. 

List  S  of  a  box  b  will  be  denoted  Vi  and  is  formed  by  all  the  children  of  the  colleagues  of  6's 
parent  that  are  well  separated  from  6. 

List  S  of  a  box  b  will  be  denoted  by  Wi.  is  empty  if  6  is  a  parent  box,  and  consists  of  all 
descendants  of  6’s  colleagues  whose  parents  are  adjacent  to  b,  but  who  are  not  adjacent  to  b 
themselves,  if  6  is  a  childless  box.  Note  that  b  is  separated  from  each  box  w  in  IVj  by  a  distance 
greater  than  or  equal  to  the  length  of  the  side  of  w. 

List  4  of  a  box  6  will  be  denoted  by  Aj  and  is  formed  by  all  boxes  c  such  that  b  e  Wg.  Note 
that  all  boxes  in  List  4  are  childless  and  larger  than  b. 

List  5  of  a  box  b  will  be  denoted  by  Yi  and  consists  of  all  boxes  that  are  well  separated  from 
6’s  parent. 

will  denote  the  p-tcrm  multipole  expansion  about  the  center  of  6  of  the  field  created  by  all 
particles  in  r(6). 

'9ft  will  denote  the  p-term  local  expansion  about  the  center  of  box  6  of  the  field  created  by  all 
particles  located  outside  T{Ut)  U7’(n  4).  is  the  result  of  the  evaluation  of  the  expansion 

'9tt  at  a  particle  r  in  r(6). 

r»  will  denote  the  local  expansion  about  the  center  of  6  of  the  field  due  to  all  particles  in  T{Vt). 
At  will  denote  the  local  expansion  about  the  center  of  6  representing  the  field  due  to  all  charges 
located  in  T{Xh). 

at{r)  will  denote  the  the  field  at  r  €  r(6)  due  to  all  particles  in  T{Ut). 

0t{r)  will  denote  the  field  at  r  €  r(6)  due  to  all  particles  in  T{}Vb), 

3.3  Informal  description  of  the  algorithm 

The  algorithm  can  be  viewed  as  a  recursive  process  of  subdividing  the  computational  cell  into 
increasingly  finer  meshes  (see  Figures  2-3).  For  a  fixed  box  6  at  level  /,  the  computational  cell 
is  partitioned  into  five  subsets,  Uf,  V*,  IVj,  A'j  and  I »,  and  the  following  procedure  is  applied 
to  the  sets  -if  particles  T{Ut),  T(Vt)^  T{\Vt),  T{Xt)  and  T'(li). 

1.  For  each  childless  box  6  we  combine  the  particles  in  r(6)  by  means  of  Theorem  2.1  to 
form  a  multipole  expansion  For  each  parent  box  B  we  use  Lemma  2.2  to  merge  the 
multipole  expansions  of  its  children  6i,  62,  63  ,6^  into  the  expansion  $£. 

2.  The  interactions  between  particles  in  r(6)  and  T{Ut)  are  computed  directly.  For  each 
particle  r  €  T{b),  the  result  of  these  calculations  is  at(r). 

3.  We  use  Lemma  2.3  to  convert  the  multipole  expansion  of  each  box  in  into  a  local 
expansion  about  the  center  of  6,  and  add  the  resulting  expansions  to  obtain 

4.  For  every  particle  r  in  6,  we  compute  the  field  /?»(r)  due  to  all  particles  in  T{\Vt)  by 
evaluating  the  p-term  multipole  expansions  of  each  box  w  in  M't  at  r,  and  adding 
them  up. 
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5.  We  convert  the  field  of  each  particle  in  T’(A'j)  into  a  local  expansion  about  the  center  of 
box  6  (see  Remark  2.1),  and  add  up  the  resulting  expansions  obtaining  Aj. 

6.  We  shift  the  center  of  the  local  expansion  Tb  of  fc’s  parent  B  to  the  centers  of  6  and  the 
other  children  of  B  by  means  of  Lemma  2.4.  We  add  the  local  expansion  obtained  to  Ft. 

7.  For  each  box  6,  we  evaluate  the  sum  of  the  local  expansions  r6  and  A^  at  every  particle 
r  in  b  and  add  the  result  to  on{r)  and  ^^{r)  obtaining  the  field  at  r. 

Remark  3.1.  Note  that  in  the  above  procedure  we  never  explicitly  evaluate  the  interactions 
between  particles  in  T{b)  and  these  in  r(l*).  Indeed,  since  all  boxes  in  Vj  are  well  separated 
from  6’s  parent,  the  interaction  between  T{Y\)  and  T{b)  have  been  accounted  for  during  steps 
3  and  5  at  a  coarser  level. 

3.4  Formal  description  of  the  algorithm 

Algorithm 

Commer.t  [Choose  main  parameters] 

Choose  precision  e  to  be  achieved.  Set  the  number  of  terms  in  all  expansions  to  p  se  Iog2(f ). 
Choose  the  maximum  number  s  of  particles  in  a  childless  box. 

Stage  1. 

Comment  [Refine  the  computational  cell  into  a  hierarchy  of  meshes.) 

do  1  =  1,2,  •  •  • 
do  6,'  €  Bi 

if  bi  contains  more  than  s  particles  then 

subdivide  6,-  into  four  boxes,  ignore  the  empty  boxes  formed, 
add  the  non-empty  boxes  formed  to  B/+i. 
end  if 
end  do 
end  do 

Comment  [We  denote  by  nlev  the  highest  level  of  refinement,  and  by  nbox  the  total  number 
of  boxes  formed  in  Stage  l.j 

Stage  2. 

Comment  [For  every  box  b  at  every  level  f ,  form  a  multipole  expansion  representing  the  field 
outside  b  due  to  all  the  particle*  contained  in  b.) 

Step  2.1 

Comhient  [For  each  childless  box  b,  use  Theorem  2.1  to  combine  all  charges  inside  b  to  obtain 
the  multipole  expansion  about  the  center  of  b.) 
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do  !=:l,nbox 

if  bi  is  a  childless  box,  use  Theorem  2.1  to  form  a  j9-term  expansion  <^4, 
represeniing  the  field  outside  6,- due  to  all  charges  located  in  6,-. 
end  do 

Step  2.2 

Comment  [  For  each  parent  box  6,  use  Lemma  2.2  to  obtain  the  multipole  expansion  ^>4 
by  shifting  the  centers  of  the  expansions  of  6’s  children  to  h’s  center,  and  adding  the  resulting 
expansions  together.] 

So  l=nlev-l,l,-l 
do  bi  €  Bi 

if  bi  is  a  parent  box  then 

use  Lemma  2.2  to  shift  the  center  of  each  of  6,’s  child  box’s  expansion  to 
6,'s  center.  Add  the  resulting  expansions  together  to  obtain  the  expansion  $4,. 
end  if 
end  do 
end  do 


Stage  3. 

Comment  (For  all  particles  in  each  childless  box  6,  compute  the  interactions  with  all  particles 
in  T{Ut,)  directly.) 

do  i=i,nbox 

if  bi  is  childless  then 

for  each  particle  r  in  6,-,  compute  the  sum  a4(r)  of  the  interactions  between  r 
and  all  particles  in  T{Ui.). 
end  if 
end  do 


Stage  4. 

Comment  [  For  each  box  6,  use  Lemma  2.3  to  convert  the  multipole  expansions  of  all  boxes 
in  V4  into  local  expansions  about  the  center  of  box  b.) 

do  i— l,nbox 
do  bj  €  V^^- 

Convert  multipole  expansion  auout  6y’s  center  into  a  local  expansion 
about  6,’’s  center  using  Lemma  2.3.  Add  the  resulting  expansions  to  obtain 
end  do 
end  do 


Stage  5. 
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Comment  [For  each  childless  box  b,  evaluate  the  multipole  expansions  of  all  boxes  in  IT  j  at 
every  particle  position  in  6.] 

do  i=l,nbox 

if  6,-  is  childless  then 

Evaluate  the  multipole  expansion  of  each  box  bj  €  iVt,  to  obtain 
0t,(r)  for  every  particle  r  in  box  b,-. 
end  if 
end  do 

Stage  6. 

Comment  [For  each  box  b,  use  Lemma  2,3  and  Remark  2.2  to  form  local  expansions  about 
the  center  b  representing  the  field  due  to  all  particles  in  r(A'6).] 

do  i=l,nbox 

Convert  the  field  of  every  particle  in  T(XtJ  into  a  local  expansion  about  the  center  of  6 
end  do 

Stage  7. 

Comment  [Use  Lemma  2.4  to  shift  the  centers  of  local  expansions  of  parent  boxes  to  the 
centers  of  their  children.] 

do  1— l,nlev-l 
do  6,'  €  Bi 

if  6,  is  a  parent  box  then 

by  using  Lemma  2.4,  shift  the  center  of  expansion  Fj,  to  the  center 
of  each  of  6,  ’s  children  bj.  Add  the  resulting  expansion  to  Fty. 
end  if 
end  do 
end  do 

Stage  8. 

Comment  [  For  each  childless  box  b,  obtain  as  the  sum  of  local  expansions  r6  and  For 
each  particle  r  in  a  childless  box  6,  evaluate  and  obtain  the  field  at  r  by  adding  'i'fc(r), 

ot(r)  and  together,] 

do  issl,nbox 

if  bi  is  childless  then 

Compute  'i'j.  =  r».  +  Aj,.. 

For  each  particle  r  in  6,-,  evaluate 

Add  and  to  obtain  the  field  at  r’s  position, 

end  if 
end  do 
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3.5  Complexity  analysis 

Stage  Operation  Explanation 

number  count 

Stage  1  Np  Each  particle  is  assigned  to  a  box  at  every  level. 

There  are  at  most  p  levels  of  refinement. 

Stage  2 

Step  2.1  Np  Each  particle  contributes  to  the  p-term  expansion  of 

one  childless  box. 

Step  2.2  ^p^N/s  The  center  of  the  expansion  of  each  box  is  shifted  to 

the  center  of  the  parent  box.  The  number  of  boxes 
is  bounded  by  5pN/s  (see  Lemma  A. 5),  and  each  shift 
requires  p*/2  work  (see  Lemma  2.2). 

Stage  3  ISp.Vs  Each  childless  box  b  contains  less  than  s  particles 

and  the  work  required  to  compute  all  interactions  between 
particles  in  two  boxes  is  s^/2  when  Newton’s  third 
law  is  used.  The  number  of  boxes  in  all  List  I's  is  bounded 
by  SGpN/s  (see  Lemmas  A.l  and  A.4). 

Step  4  SOp^N/s  For  each  box.  List  2  has  no  more  than  32  entries  (Lemma  A.2). 

There  are  at  most  SpA’/s  boxes  (Lemma  A. 5)  and 
each  shift  requires  p*/2  work  (Lemma  2.3). 

Step  5  32p^.V  Each  childless  box  b  contains  less  than  s  particles. 

The  interactions  of  all  particles  in  b  and  a 

box  in  IV4  require  ps  work.  The  total  number  of 

boxes  in  List  3  is  bounded  by  32pAys  (Lemma  A.3  and  A.4). 

Step  6  32p^A’'  Each  box  in  A'j  contains  less  than  s  particles. 

The  interactions  between  all  particles  in  a  box  in  A'» 

and  box  6  require  ps  work.  The  total  number  of 

boxes  in  List  4  is  bounded  by  Z2pN/s  (Lemma  A.3  and  A.4). 

Stage  7  lOp^N/a  Each  box  has  at  most  four  children. 

There  are  less  than  bpN/s  boxes  (Lemma  A.5) 
and  a  shift  requires  p^/2  work  (Lemma  2.4). 

Stage  8  Np+  N  A  p-term  expansion  is  evaluated  at  each  particle  position. 

The  sums  require  an  extra  N  work. 
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Summing  up  the  CPU  times  for  all  stages  above,  we  obtain  the  following  time  estimate: 

T  =  N  ■  {92.bap^/s  +  64bp^  +  IScps  +  3dp+ e),  '21) 


where  the  coefficients  a,b,c^d^e  depend  on  the  computer  system,  language,  implementation, 
etc.  However,  the  parameter  a  (maximum  permitted  number  of  particles  in  a  childless  box) 
in  (21)  is  not  determined  by  the  problem  and  can  be  choosen  so  as  to  minimize  the  resulting 
CPU  time  estimate.  Differentiating  (21)  with  respect  to  s,  we  obtain: 


®»nin 


m  =  '/• 


92.5a 

18c 


•P 


and 


Tmin  =  A'  •  (op^  +  I3p  +  ‘j)  =  N  •  (alog^(f)  +  0\  log2(f)|  +  Tf), 


(22) 

(23) 


with  the  constants  o,/3,'y  determined  by  the  computer  system,  language,  implemantation,  etc. 

The  storage  requirements  of  the  algorithm  are  determined  by  the  number  of  non-empty 
boxes  which  is  bounded  by  5pN/s.  For  each  box  we  store  the  coefficients  of  a  p-term  multipole 
expansion  and  a  p-term  local  expansion.  The  positions  and  charges  of  each  particle  also  have 
to  be  stored.  Therefore  the  storage  requirements  are  of  the  form: 


S  =  iN'  •  (10/p/s  -f  3i7)  =  .V  •  (10/1  log2(0l/«  +  3p),  (24) 

where  the  coefficients  f,g  depend  on  the  computer  system,  language,  implementation,  etc. 

4.  Numerical  results 

A  computer  program  using  the  algorithm  described  in  the  preceding  section  has  been  imple¬ 
mented,  and  numerical  experiments  have  been  performed  on  a  VAX-8600. 

To  evaluate  the  robustness  of  the  adaptive  scheme,  we  considered  a  variety  of  particle 
distributions.  For  each  distribution,  the  corresponding  fields  were  computed  in  four  ways:  by 
the  algorithm  of  the  present  paper,  by  the  algorithm  described  in  [7],  and  directly  in  single  and 
double  precision.  The  direct  calculation  of  the  field  in  double  precision  was  used  as  a  standard 
for  comparing  the  relative  accuracies  of  the  other  three  methods.  In  these  experiments,  the 
number  of  particles  varied  between  100  and  25600,  with  charge  strengths  randomly  assigned 
between  zero  and  one. 

The  results  are  summarized  in  Tables  1, 2,  3,  and  4.  The  first  column  of  each  table  contains 
the  number  of  particles  N  for  which  calculations  have  been  performed.  In  the  remaining 
columns,  the  upper  case  letters  T,  E  and  S  are  used  to  denote  the  corresponding  computational 
time,  error  and  storage,  with  the  subscripts  alg,  uni  and  dir  referring  to  the  adaptive  algorithm, 
the  algorithm  described  in  [7],  and  the  direct  (single-precision)  calculation  respectively.  More 
specifically,  columns  2  through  4  show  the  times,  in  seconds,  required  to  compute  the  field  by 
the  three  methods.  The  errors  JS!*!,,  Euni  and  Ejir  for  the  adaptive,  non-adaptive  and  direct 
method,  respectively,  are  presented  in  the  next  three  columns.  They  are  defined  by  the  formula 
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1/2 


[  EJL.i/.f  j 

where  /.-.is  the  value  of  the  field  at  the  i-th  particle  position  obtained  by  direct  calculation  in 
double  precbion  and  /,•  b  the  result  obtained  by  one  of  the  three  methods  being  studied.  The 
last  two  columns  of  the  tables  contain  the  storage  requirements  Saig  and  Suni ,  in  single-precision 
words,  for  the  two  methods  based  on  multipole  expansions. 

Remark  4.1.  For  the  tests  involving  12800  and  25600  particles,  it  was  not  considered  practical 
to  use  the  direct  method  to  calculate  the  fields  at  all  particle  positions,  since  this  would  require 
prohibitive  amounts  of  CPU  time  without  providing  much  useful  Information.  Therefore,  we 
have  performed  the  direct  calculations  in  double  precision  for  only  100  of  the  particles,  and 
used  these  results  to  evaluate  the  relative  accuracies.  The  corresponding  values  of  Tdir  were 
estimated  by  extrapolation. 

For  the  first  set  of  tests,  the  positions  of  the  charged  particles  w'ere  randomly  distributed 
in  a  square,  and  the  resulting  particle  density  was  roughly  uniform  (Figure  6),  The  number  of 
terms  in  the  expansions  was  set  to  20,  and  the  maximum  number  of  particles  in  a  childless  box 
was  set  to  30. 

In  the  second  set  of  experiments,  the  charged  particles  were  distributed  along  a  curve 
(Figure  7).  The  number  of  terms  in  the  expansions  was  set  to  17  and  the  maximum  number  of 
particles  in  a  childless  box  was  set  to  30. 

The  third  set  of  numerical  experiments  was  performed  on  extremely  non-uniform  distribu¬ 
tions  of  particles  (Figures  8).  A  fifth  of  the  N  particles  were  randomly  assigned  in  a  square  of 
area  one.  Two  fifths  were  randomly  distributed  about  the  center  of  the  square  in  a  circle  of 
radius  0.003.  The  rest  of  the  particles  were  assigned  positions  inside  a  circle  of  radius  0.5  with 
a  density  inversely  proportional  to  the  square  of  the  distance  from  the  center.  The  number  of 
terms  in  the  expansions  was  set  to  17  and  the  maximum  number  of  particles  per  childless  box 
was  set  to  30. 

In  the  last  set  of  experiments,  half  of  the  particles  were  distributed  along  a  curve  similar 
to  that  of  the  second  set  of  experiments  and  the  rest  of  the  particles  were  distributed  inside 
four  circles  w’ith  a  density  inversly  proportional  to  the  square  of  the  distance  to  the  centers  of 
the  circles  (Figure  9).  The  number  of  terms  in  the  expansions  w'as  set  to  17  and  the  maximum 
number  of  particles  per  childless  box  was  set  to  30. 

The  following  observations  can  be  made  from  Tables  1,  2,  3  and  4,  where  the  results  of  the 
experiments  described  above  arc  summarized. 

1.  The  accuracies  of  the  results  obtained  by  the  algorithms  using  multipole  expansions  are 
in  agreement  with  the  error  bounds  given  in  (7), (14)  and  (18).  For  the  most  part,  the 
fast  methods  are  slightly  more  accurate  than  the  direct  calculation. 

2.  In  all  cases,  the  actual  CPU  time  requirements  of  the  adaptive  algorithm  grow  linearly 
with  N.  The  CPU  time  requirements  of  the  non-adaptive  counterpart  grow’  linearly, 
except  for  extremely  non-uniform  distributions  (see  Tables  3-4). 
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3.  Even  for  uniform  distributions  of  charges,  the  adaptive  aigoritiim  is  about  30%  faster 
than  the  non-adaptive  one. 

4.  The  storage  requirements  of  both  fast  algorithms  are  roughly  proportional  to  the  num¬ 
ber  of  particles  involved  in  the  simulations.  The  storage  requirements  of  the  adaptive 
algorithm  are  about  four  times  less  than  those  of  the  non-adaptive  version. 

5.  By  the  time  the  number  of  particles  reaches  25600,  the  adaptive  algorithm  is  about  100 
times  faster  than  the  direct  version  for  the  case  of  the  uniform  distribution  (see  Table  1). 
When  the  charges  are  situated  on  a  curve,  the  adaptive  scheme  is  roughly  200  times  faster 
than  the  direct  algorithm,  and  about  3  times  faster  than  the  fast  non-adaptive  scheme(see 
Table  2).  For  the  highly  non-uniform  case  (see  Table  3),  the  adaptive  algorithm  is  slightly 
more  efficient  than  for  the  uniform  distribution.  The  non-adaptive  scheme  displays  an 
almost  quadratic  growth  of  CPU  time  with  A',  and  is  about  25  times  slower  than  its 
adaptive  counterpart  by  the  time  X  =  25600. 

6.  Even  for  as  few  as  1600  particles,  the  adaptive  algorithm  is  about  ten  times  faster  than 
the  direct  calculation. 

7.  The  performance  of  the  algorithm  does  not  depend  on  the  shape  of  the  region  where  the 
charges  are  distributed  (see  Table  4.) 


5.  Conclusions 

/ 

An  adaptive  algorithm  ha.s  been  constucted  for  the  rapid  evaluation  of  the  potentials  and  force 
fields  due  to  large4cale  ensembles  of  particles  of  the  type  encountered  in  plasma  physics,  molec¬ 
ular  dynamics,  fluid  dynamics  (the  vortex  method),  and  celestial  mechanics.  The  algorithm  is 
applicable  whenever  the  fields  to  be  evaluated  are  Coulombic  or  gravitational  in  nature,  and 
yields  the  potentials  to  within  round-off  error. 

The  asymptotic  CPU  time  estimate  for  the  algorithm  is  of  the  order  0(,Vj,  where  N  is  the 
number  of  particles  in  the  simulation,  and  this  e.jtimate  is  independent  of  the  statistics  of  the 
charge  distribution.  Our  numerical  experiments  indicate  a  tendency  of  the  scheme  to  be  more 
efficient  for  non-uniform  distributions  than  for  uniform  ones.  The  storage  requirements  of  the 
V  algorithm  are  of  the  order  0{N),  do  not  depend  on  the  statistics  of  the  distribution,  and  tend 
I  to  be  quite  acceptable  even  for  very  large  numbers  of  particles. 

'“"•^  'Tfrthr  prw nt  paper,  a  two-dimensional  version  of  the  algorithm  is  described.  Generalizing 
it  to  the  three-dimensional  case  b  fairly  straightforward,  and  will  be  reported  at  a  later  date. 

The  authors  would  like  to  thank  Professor  M.  H.  Schultz  for  several  useful  discussions  and 
for  hb  interest  and  support. 
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Appendix 

In  this  Appendix,  we  prove  several  combinatorial  lemmas  that  are  used  in  Subsection  3,5 
to  estimate  the  complexity  of  the  algorithm  described  in  Subsections  3.3  and  3.4.  We  begin  by 
introducing  some  additional  notation. 

Given  a  subdivision  5  of  the  computational  cell  and  a  childless  box  b  in  S,  we  will  denote 
by  Sb  the  subdivision  obtained  from  S  by  subdividing  h  into  4  equal  boxes,  and  refer  to  the 
process  of  obtaining  from  5  as  an  elementary  refinement  of  5. 

For  a  subdivision  S  of  the  computational  cell,  we  will  denote  by  Bs  the  set  of  all  boxes  in  5. 
Cs  will  denote  the  subset  of  Bs  consisting  of  all  childless  boxes,  i.e.  boxes  that  are  non-empty 
and  not  subdivided. 

Fs  will  denote  the  subset  of  Bs  consisting  of  all  non-empty  boxes. 

Ds  is  the  subset  of  Bs  consisting  of  all  empty  boxes  that  have  a  childless  brother. 

Gs  —  Cs  U  Ds  is  the  subset  of  Bs  consisting  of  all  boxes  b  such  that  b  is  either  childless  or  has 
a  childless  brother. 

For  any  set  of  boxes  A,  X (A)  will  denote  the  number  of  boxes  in  A. 

Lemma  A.l:  For  any  subdivbion  S  of  the  computational  cell 

T,  (25) 

»€C5 

Proof:  We  will  prove  the  lemma  by  combining  the  following  three  observations; 

a)  Inequality  (25)  holds  for  the  undivided  computational  cell. 

b)  Any  subdivision  of  the  computational  cell  can  be  obtained  by  a  sequence  of  elementary 
refinements  of  the  computational  cell. 

c)  If  an  elementary  refinement  b  applied  to  a  subdivision  for  which  (25)  holds,  it  also  holds  in 
the  refined  subdivbion. 

The  statements  a)  and  b)  above  are  obvious,  and  following  b  a  proof  of  c). 

Consider  a  subdivbion  S  of  the  computational  cell  such  that  (25)  is  true  for  5,  and  a  box 
6  such  that  6  €  Cs.  Obviously 

A'(CsJ  =  A^(Cs))  +  3,  (26) 

and  we  will  denote  by  I/*  and  the  Lbt  I’s  of  b  with  respect  to  5  and  5*,  respectively.  The 
following  obvious  observations  can  be  made  about  the  List  I’s  of  b  and  its  children: 

1)  For  any  box  c  €  Cs,  if  6  6  t/e  then  c^U^. 

2)  Each  child  of  6  has  its  three  brothers  in  its  Lbt  1. 

3)  In  the  subdivbion  5»,  b  b  not  childless  and  Ul  b  empty. 

4)  Each  box  c  in  b  in  the  List  1  of  at  least  one  child  of  6. 

5)  The  number  of  boxes  of  that  are  in  the  Lbt  Is  of  two  children  of  6  b  bounded  by  8. 

It  immediately  follows  from  observations  l)-5)  above  that 

E  E  =  ®  +  *  +  +  (27) 

p€Cs^  f€Cs 
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(28) 


and  combining  (26)  and  (27)  we  obtain 

L«mma  A. 2:  For  any  subdivbion  5  of  the  computational  cell 

E  5  32  •  tflfs).  (29) 

iefs 

Proof:  Consider  a  arbitrar}*  subdivision  S  of  the  computational  cell,  a  box  ce  Fs  and  its  parent 
box  6.  is  a  subset  of  the  children  of  6*s  colleagues,  the  maximum  number  of  colleagues  of  c 
(or  any  other  box)  b  eight,  and  each  colleague  can  have  four  children.  Therefore,  the  number 
of  elements  in  b  bounded  by  32. 

Lemma  A.3:  For  any  subdivbion  5  of  the  computational  cell 

E  *''’(”*)  =  E  ^  S  •  (30) 

€€Cs  ^Fs 

Proof:  The  first  part  of  the  Lemma  is  a  direct  consequence  of  the  definition  of  List  4  (see 
Subsection  3.2):  If  a  box  b  belongs  IV,,  then  e  belongs  to  .Y».  Now,  consider  an  arbitrary  box 
c  €  Fs,  and  its  parent  box  b.  The  number  of  colleagues  of  b  b  certainly  bounded  by  8.  We 
will  denote  by  2*  the  set  of  all  childless  boxes  which  are  adjacent  to  b  and  whose  size  is  greater 
than  or  equal  to  that  of  b.  The  number  of  boxes  in  2*  is  bounded  by  8,  since  each  box  in  2^ 
contains  at  least  one  of  the  eight  colleagues  of  6,  and  no  two  such  boxes  can  contain  the  same 
colleague.  The  second  part  of  the  lemma  now  follows  from  the  the  obvious  observation  that 
\Vc  C  2*. 

A.4:  For  any  subdivbion  5  of  the  computational  cell  produced  by  the  algorithm  of 
Section  3, 

A'(Cs)<N’((?s)<4  p— .  (31) 

Proof:  Each  parent  box  b  at  lev'el  /  contains  more  than  s  particles  (otherwise,  it  would  not 
have  been  subdivided  any  further).  Therefore,  the  total  number  of  parent  boxes  at  level  /  is 
bounded  by  N/s.  Each  of  these  boxes  can  not  have  more  than  4  children,  and  consequently  the 
number  of  boxes  in  Gs  at  any  level  /  b  bounded  by  4iV/«.  Now,  the  conclusion  of  the  lemma 
follows  from  Observation  3.1  and  the  obvious  fact  that  S{Cs)  <  N{Gs)- 

Lemma  A.5:  For  any  subdivbion  5  of  the  computational  Cv*ll, 

N{Fs)<5’P^.  (32) 

Proof;  The  number  of  parent  boxes  at  any  level  I  b  bounded  by  N/s,  end  each  of  them  can  not 
have  more  than  4  childless  boxes  at  level  /+ 1.  Therefore,  the  sum  of  the  numbers  of  non*empty 
boxes  (all  childless  and  parent  boxes)  at  all  leveb  is  bounded  by  p  •  {N/$  +  iN/s). 
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N 

MJBM 

Tdir 

\mmm 

WMsmm 

Edir 

100 

0.15 

0.47 

0.15 

1.7  10"® 

■TPlTtBl 

1.7  10"® 

866 

4179 

200 

0.43 

0.65 

0.61 

9.3  10"^ 

dEB 

4.4  10“^ 

2503 

5479 

400 

1.01 

1.94 

2.47 

7.0  10“* 

6.4  10"’ 

3763 

16847 

800 

2.45 

2.78 

10.27 

4.1  10"^ 

4.0  10"^ 

4.7  10"’ 

11203 

22047 

1600 

5.37 

8.56 

42.35 

3.7  10"^ 

4.2  10"^ 

5.4  10"’ 

15923 

67519 

3200 

10.60 

11.80 

152.95 

5.0  10"^ 

5.3  10"'^ 

8.7  10"’ 

44423 

88319 

6400 

23.38 

33.49 

601.18 

7.0  10“' 

5.4  10"’ 

1.3  10“® 

65907 

270207 

12800 

45.34 

48.02 

2433.20 

6.0  10"^ 

4.9  10"’ 

1.6  10"« 

176631 

353407 

25600 

96.72 

137.68 

9694.45 

8.3  10"* 

8.9  10"’ 

2.2  10"« 

268723 

1080959 

Table  1:  Uniformly  distributed  particles, 
p  =  20  and  5  =  30 


N 

mm 

mm 

“Edir 

Ealf 

Edit 

IHili 

^EIBi 

100 

0.11 

0.38 

0.16 

IffflTtti 

l(Kh 

1149 

2C0 

0.30 

0.54 

0.57 

nffig 

2694 

0.64 

1.31 

2.29 

5.6  10"® 

BPS  [Bfl 

5103 

15827 

800 

1.46 

3.13 

9.30 

9.4  10“* 

9.5  10"® 

9.5  10"® 

10133 

21027 

1600 

2.6G 

5.94 

37.41 

2.0  10"* 

2.0  10"® 

2.0  10"® 

19241 

63427 

3200 

5.93 

12.50 

149.21 

7.8  10"® 

8.7  10"® 

8.8  10"® 

40055 

84227 

’6400 

12.42 

29.66 

597.95 

4.2  10"® 

4.2  10"® 

4.2  10"® 

84429 

12800 

25.11 

79.47 

2425.48 

8.7  10"® 

8.7  10"® 

t.8  10"® 

167421 

■?  SwSB 

25600 

47.53 

152.07 

9581.20 

8.9  10"® 

9.1  10"® 

8.9  10"® 

332927 

1015427 

l^ble  2;  Particles  distributed  on  a  curve. 
p=  17  and  s  =  30 


N 

\mm 

flE2E9l 

mmm 

Ear 

•S'uni 

100 

0.19 

0.45 

2.7  lO"*" 

mwmaM 

2508 

3927 

200 

0.48 

0.74 

6.9  10-* 

7.6  10-® 

4014 

5227 

400 

1.13 

2.26 

1.9  10"* 

9.0  10-® 

8307 

15827 

800 

2.25 

5.15 

4.3  10-« 

3.7  10-® 

13353 

21027 

1600 

5.09 

16.17 

37.74 

2.4  !0-® 

2.1  10-® 

25588 

63427 

3200 

9.98 

50.23 

149.86 

3.7  10-* 

1.4  10-® 

1.7  10-® 

46806 

84227 

6400 

21.80 

177.13 

606.14 

5.8  10-® 

4.0  10-® 

5.9  10”® 

90505 

253827 

12800 

41.93 

663.21 

2420.33 

4.0  10-* 

4.2  10”® 

186226 

337027 

25600 

90.05 

2317.93 

9622,63 

2.9  10-® 

3.0  10"® 

4.0  10”® 

373639 

1015427 

Table  3:  Highly  non-uniform  distribution  of  particles, 
p  =  17  and  «  =  30 


.v 

MUM 

mHH 

'Edir 

wEnm 

WEBStk 

Edir 

11091 

•^uni 

0.15 

0.43 

0.15 

4.3  10”^ 

EQmBil 

1145 

0.39 

0.68 

0.59 

3.3  10”® 

3224 

0.84 

1.69 

2.31 

RfifB 

7.1  10”® 

8.1  10”® 

6939 

15827 

800 

2.11 

5.03 

9.39 

4.3  10”® 

4.3  10”® 

4.3  10”® 

13406 

21027 

1600 

4.35 

11.34 

37.74 

9.2  10”® 

9.2  10”® 

9.2  10”® 

24913 

63427 

3200 

9.16 

30.85 

153.76 

1.1  10”® 

1.1  10”® 

1.1  10”® 

48902 

84227 

6400 

19.22 

48.62 

611.82 

5.4  10”® 

5.5  10”® 

5.4  10”® 

96153 

253827 

12800 

37.92 

155.75 

2440.90 

2.1  10”® 

2.0  10”® 

2.1  10”® 

194377 

337027 

25600 

80.02 

248.90 

' 

9798.34 

4.4  10”® 

4.4  10”® 

4,5  10”® 

388624 

1015427 

Table  4:  Non-uniform  dbtribution  of  particles  in  a  region  of  complicated  shape, 

p  =  17  and  a  =  30 


Figure  6:  25600  uniformly  located  charges  in  the  computational  cell. 


Figure  8:  Righlj  non-uniform  distribution  of  25600  chsrges. 


Figure  9:  A  iieii>itBifonD  dittribatSeo  of  25600  particles  in  a  region  of  complicated  shape. 


